Load packages

## add 'developer' to packages to be installed from github
x <- c("data.table", "lubridate", "devtools", github = "maRce10/warbleR", "readxl", "ranger", "caret", "e1071", "pbapply", "viridis", "ggplot2", "DT", "kableExtra", "rlang", "Sim.DiffProc", "soundgen", "brms", "ggrepel", "posterior", "ggridges", "tidybayes", "purrr", bioconductor = "multtest", "metap"#, "markovchain", "igraph", "TraMineR", "spgs"
       )

source("~/Dropbox/R_package_testing/sketchy/R/load_packages.R")
load_packages(x, quite = TRUE)

Functions and global parameters

warbleR_options(wl = 300, parallel = 1, bp = "frange", fast = TRUE, threshold = 15, ovlp = 20)

opts_knit$set(root.dir = "..")


theme_set(theme_classic(base_size = 24))

# set evaluation false
opts_chunk$set( fig.width = 7, fig.height = 4, warning = FALSE, message = FALSE, tidy = FALSE)

cols <- viridis(10, alpha = 0.7)

source("~/Dropbox/R_package_testing/brmsish/R/html_summary.R")
source("~/Dropbox/R_package_testing/brmsish/R/contrasts.R")
source("~/Dropbox/R_package_testing/brmsish/R/helpers.R")

# model parameters
chains <- 4
iter <- 10000

prior <- c(prior(normal(0, 10), "b"), prior(normal(0, 50), "Intercept"),
    prior(student_t(3, 0, 20), "sd"), prior(student_t(3, 0, 20), "sigma"))

Read data

clls <- readRDS("./data/processed/curated_extended_selection_table_inquiry_calls_2020_&_2021.RDS")

sub_clls <- clls[!duplicated(clls$org.sound.files), ]
sub_clls$file.duration <- sub_clls$sound.file.samples / (sub_clls$sample.rate * 1000)

metadat <- read.csv("./data/processed/metadata_inquiry_calls_2020_&_2021.csv", stringsAsFactors = FALSE)

caps <- as.data.frame(read_excel("./data/raw/Proyecto MPI enero 2020_2.xlsx", sheet = "Capturas"))

# read acoustic parameter data
acous_param_l <- readRDS("./data/processed/acoustic_parameters_all_groups_specific_warbler_acoustic_measurements_curated_data_2020_&_2021.RDS")

# read as RDS
agg_pred <- readRDS("./data/processed/predicted_individual_in_group_flights_2020_&_2021.RDS")

# read diagnostics
diagnostics <- readRDS("./data/processed/random_forests_diagnostics_solo_flight.RDS")
# remove groups with individuals with no solo flight

diagnostics$all_with_solo <- sapply(1:nrow(diagnostics), function(x) 
{
    Y <- acous_param_l[[diagnostics$group[x]]]
  Y$sound.files <- sapply(Y$sound.files, function(x) strsplit(x, ".wav")[[1]][1])
  Y$year.audio <- sapply(Y$sound.files, function(x) sub_clls$year.audio[gsub(".wav", "", sub_clls$org.sound.files) == x][1])
    
  indivs <- strsplit(na.omit(Y$indiv[Y$year.audio == diagnostics$group_flight_files[x]])[1], split = "\\|")[[1]]

  call_per_indiv <- sapply(indivs, function(y) sum(clls$Individuo == y & clls$Experimento == "vuelo solo"))
  
  if (any(call_per_indiv == 0)) return("missing solo")  else 
    return("OK")
  }
)

diagnostics <- diagnostics[diagnostics$all_with_solo == "OK", ]

agg_pred <-agg_pred[agg_pred$group %in% diagnostics$group, ]

# get summary by group
summary_by_group_list <- lapply(unique(agg_pred$group), function(x) {
    # print(x)
  Y <- agg_pred[agg_pred$group == x, ]

#total number of calls above lowest group threshold for true positives
  above_threshold_calls <- sum(Y$max_prob > diagnostics$min_prob_threshold[diagnostics$group == x])

#proportion of calls
  prop_above_calls <- above_threshold_calls / nrow(Y)

  return(data.frame(group = x, experiment = diagnostics$experiment[diagnostics$group == x], total_calls = nrow(Y), above_threshold_calls = above_threshold_calls, prop_above_calls, n_individuals = diagnostics$n_indiv[diagnostics$group == x]))
  
})

# calls per individual and group
calls_by_indiv_group <- aggregate(sound.files ~ group + pred_indiv, agg_pred, length)

names(calls_by_indiv_group) <- c("Group", "Individual_ID", "Number_of_calls")

calls_by_indiv_group <- calls_by_indiv_group[order(calls_by_indiv_group$Group), ]

write.csv(calls_by_indiv_group, "./data/processed/calls_by_indiv_group.csv", row.names = FALSE)

summary_by_group <- do.call(rbind, summary_by_group_list)

groups_by_cat <- aggregate(group ~ experiment, summary_by_group, FUN = function(x) length(unique(x)))

#print pretty table
df3 <- knitr::kable(groups_by_cat, row.names = FALSE, escape = FALSE, format = "html", digits = 2)

kable_styling(df3, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, font_size = 18)
experiment group
mixed 33
regular 26

Group count by number of individuals and experiment

groups_by_cat_n <- as.data.frame.matrix(table(summary_by_group$experiment, summary_by_group$n_individuals))

#print pretty table
df4 <- knitr::kable(groups_by_cat_n, row.names = TRUE, escape = FALSE, format = "html", digits = 2)

kable_styling(df4, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, font_size = 18)
2 3 4 5 6 7
mixed 12 10 7 4 0 0
regular 2 7 8 5 2 2

Call rate

# get group and solo call counts per individual, call rate in call/min
call_rate_list <- pblapply(1:nrow(diagnostics), function(x){

  Y <- acous_param_l[[diagnostics$group[x]]]
  Y$sound.files <- sapply(Y$sound.files, function(x) strsplit(x, ".wav")[[1]][1])
  Y$year.audio <- sapply(Y$sound.files, function(x) sub_clls$year.audio[gsub(".wav", "", sub_clls$org.sound.files) == x][1])
  
  indivs_in_group <- strsplit(Y$indiv[Y$year.audio == diagnostics$group_flight_files[x]], split = "\\|")[[1]] 
  
  calls_per_sound_file_list <- lapply(na.omit(c(indivs_in_group, "group")), function(y){ 
    
      if (y != "group"){
      tab <- table(Y$year.audio[Y$indiv == y])
      if (length(tab) == 0) {tab <- 0
      names(tab) <- y}
      } else{
        tab <- sum(Y$year.audio == diagnostics$group_flight_files[x], na.rm = TRUE)
    names(tab) <- diagnostics$group_flight_files[x]
        }
    
    df <- data.frame(group = diagnostics$group[x], experiment = diagnostics$experiment[x], year.audio = names(tab), indiv = y, calls = as.vector(tab), flight.time = sapply(names(tab), USE.NAMES = FALSE, function(x) (round(as.numeric(metadat$Tiempo.de.vuelo[metadat$year.audio == x][1]) * 60 / 0.04166667))), file.duration = sapply(names(tab), function(x) (sub_clls$file.duration[sub_clls$year.audio == x][1])))
    
    return(df)
    }
  )
  
  calls_per_sound_file <- do.call(rbind, calls_per_sound_file_list)
  
  calls_per_sound_file$rate <- calls_per_sound_file$calls / calls_per_sound_file$file.duration * 60 # calls 

  return(calls_per_sound_file)
    })

call_rate_by_group <- do.call(rbind, call_rate_list)

rownames(call_rate_by_group) <- 1:nrow(call_rate_by_group)

#  add call rate to those missing
call_rate_by_group$rate <- ifelse(!is.na(call_rate_by_group$rate), call_rate_by_group$rate, call_rate_by_group$calls / call_rate_by_group$flight.time)

# add sex
call_rate_by_group$sex <- sapply(call_rate_by_group$indiv, function(x) if(x != "group") na.exclude(caps$Sexo[caps$Murci == x])[1] else NA, USE.NAMES = FALSE)

call_rate_by_group$sex <- ifelse(call_rate_by_group$sex == "m", "Male", "Female")

# add age
call_rate_by_group$age <- sapply(call_rate_by_group$indiv, function(x) if(x != "group") na.exclude(caps$Edad[caps$Murci == x])[1] else NA, USE.NAMES = FALSE)

call_rate_by_group$age <- ifelse(call_rate_by_group$age == "sa", "Sub-adult", "Adult")

# reproductive stage
call_rate_by_group$reprod.stg <- sapply(call_rate_by_group$indiv, function(x) if(x != "group") na.exclude(caps$`Estado reproductivo`[caps$Murci == x])[1] else NA, USE.NAMES = FALSE)

call_rate_by_group$reprod.stg[call_rate_by_group$reprod.stg == "p?"] <- "Pregnant"
call_rate_by_group$reprod.stg[call_rate_by_group$reprod.stg == "ne"] <- "Inactive"

call_rate_by_group$group.size <-sapply(call_rate_by_group$group, function(x) diagnostics$n_indiv[diagnostics$group == x])

# order levels
call_rate_by_group$experiment <- factor(call_rate_by_group$experiment, levels = c("regular", "mixed"))

call_rate_by_group$type <- factor(ifelse(call_rate_by_group$indiv == "group", "group", "solo"), levels = c("solo", "group"))

saveRDS(call_rate_by_group, "./data/processed/call_rate_by_group.RDS")
# get group and solo call counts per individual, call rate in call/min
call_rate_indiv_list <- pblapply(1:nrow(diagnostics), function(x){

  # print(x)
  Y <- acous_param_l[[diagnostics$group[x]]]
  
  Y$pred.id <- NA
  # add id for group flights
  Y$pred.id[Y$experiment != "vuelo solo"] <- sapply(Y$sound.files[Y$experiment != "vuelo solo"], function(x) agg_pred$pred_indiv[agg_pred$sound.files == x])
  
  Y$sound.files <- sapply(Y$sound.files, function(x) strsplit(x, ".wav")[[1]][1])
  Y$year.audio <- sapply(Y$sound.files, function(x) sub_clls$year.audio[gsub(".wav", "", sub_clls$org.sound.files) == x][1])
  
  indivs_in_group <- na.omit(strsplit(Y$indiv[Y$year.audio == diagnostics$group_flight_files[x]], split = "\\|")[[1]]) 
  
  rate_indiv_list <- lapply(indivs_in_group, function(y){ 
    # print(y)
    year.audio.group <- na.omit(unique(Y$year.audio[Y$experiment != "vuelo solo"]))
    year.audio.solo <- na.omit(unique(Y$year.audio[Y$indiv == y & Y$experiment == "vuelo solo"]))
    
    flight.time.group <-  round(as.numeric(metadat$Tiempo.de.vuelo[metadat$year.audio == year.audio.group][1]) * 60 / 0.04166667)
  
    flight.time.solo <- sapply(year.audio.solo, function(w) round(as.numeric(metadat$Tiempo.de.vuelo[metadat$year.audio == w][1]) * 60 / 0.04166667)
)  
    file.duration.group <- sub_clls$file.duration[sub_clls$year.audio == year.audio.group][1]
    file.duration.solo <- sapply(year.audio.solo, function(w) sub_clls$file.duration[sub_clls$year.audio == w][1])
    
    calls.solo <- sapply(year.audio.solo, function(w) sum(Y$year.audio == w))
    calls.group <- sum(Y$year.audio == year.audio.group & Y$pred.id == y)
    
    df <- data.frame(indiv = y, group = diagnostics$group[x], size = length(indivs_in_group), experiment = diagnostics$experiment[x], type = c("group", "solo"), year.audio = c(year.audio.group, paste(year.audio.solo, collapse = "/")), calls = c(calls.group, sum(calls.solo)), flight.time = c(flight.time.group, sum(flight.time.solo)), file.duration = c(file.duration.group, sum(file.duration.solo)))
    
    return(df)
    }
  )
  
  rate_indiv <- do.call(rbind, rate_indiv_list)
  
  rate_indiv$rate <- rate_indiv$calls / rate_indiv$file.duration * 60 # calls 

  return(rate_indiv)
    })

call_rate_indiv <- do.call(rbind, call_rate_indiv_list)

rownames(call_rate_indiv) <- 1:nrow(call_rate_indiv)

#  add call rate to those missing
call_rate_indiv$rate <- ifelse(!is.na(call_rate_indiv$rate), call_rate_indiv$rate, call_rate_indiv$calls / call_rate_indiv$flight.time)

# add sex
call_rate_indiv$sex <- sapply(call_rate_indiv$indiv, function(x) if(x != "group") na.exclude(caps$Sexo[caps$Murci == x])[1] else NA, USE.NAMES = FALSE)

call_rate_indiv$sex <- ifelse(call_rate_indiv$sex == "m", "Male", "Female")

# add age
call_rate_indiv$age <- sapply(call_rate_indiv$indiv, function(x) if(x != "group") na.exclude(caps$Edad[caps$Murci == x])[1] else NA, USE.NAMES = FALSE)

call_rate_indiv$age <- ifelse(call_rate_indiv$age == "sa", "Sub-adult", "Adult")

# reproductive stage
call_rate_indiv$reprod.stg <- sapply(call_rate_indiv$indiv, function(x) if(x != "group") na.exclude(caps$`Estado reproductivo`[caps$Murci == x])[1] else NA, USE.NAMES = FALSE)

call_rate_indiv$reprod.stg[call_rate_indiv$reprod.stg == "p?"] <- "Pregnant"
call_rate_indiv$reprod.stg[call_rate_indiv$reprod.stg == "ne"] <- "Inactive"

call_rate_indiv$group.size <-sapply(call_rate_indiv$group, function(x) diagnostics$n_indiv[diagnostics$group == x])

# order levels
call_rate_indiv$experiment <- factor(call_rate_indiv$experiment, levels = c("regular", "mixed"))

call_rate_indiv$type <- factor(ifelse(call_rate_indiv$type == "group", "group", "solo"), levels = c("solo", "group"))

saveRDS(call_rate_indiv, "./data/processed/call_rate_by_individual.RDS")
call_rate_by_group <- readRDS("./data/processed/call_rate_by_group.RDS")
call_rate_indiv <- readRDS("./data/processed/call_rate_by_individual.RDS")
call_rate_by_group <- call_rate_by_group[call_rate_by_group$type != "solo", ]
call_rate_by_group$type <- "overall.group"

call_rate_indiv$type <- as.character(call_rate_indiv$type)
call_rate_indiv$type[call_rate_indiv$type == "group"] <- "indiv.group" 

call_rate <- rbind(call_rate_by_group[, intersect(names(call_rate_by_group), names(call_rate_indiv))
], call_rate_indiv[, intersect(names(call_rate_by_group), names(call_rate_indiv))
])

# aggregate for plot
agg_rate <- aggregate(rate ~ experiment + type + group.size, data = call_rate, FUN = mean)

agg_rate$sd <- aggregate(rate ~ experiment + type + group.size, data = call_rate, FUN = sd)$rate

agg_rate$type <- factor(agg_rate$type, levels = c("solo", "indiv.group", "overall.group"))

# ggplot(agg_rate[agg_rate$group.size < 6, ], aes(x = type, y = rate)) + 
#   geom_point(color = viridis(10)[3], size = 2) +
#     geom_errorbar(aes(ymin = rate - sd, ymax = rate + sd), color = viridis(10)[3], width=.1, lwd = 1.1) +
#    geom_hline(yintercept = mean(call_rate$rate[call_rate$type == "solo"], na.rm = TRUE), lty = 2, alpha = 0.5) + 
#   labs(x = "Flight", y = "Call rate (calls / min)") + 
#   facet_grid(group.size ~ experiment) + 
#     theme(axis.text.x = element_text(angle = 45, hjust = 1))
# aggregate for plot
agg_rate_indiv <- aggregate(rate ~ experiment + type + group.size + indiv, data = call_rate_indiv, FUN = mean)

agg_rate_indiv$sd <- aggregate(rate ~ experiment + type + group.size + indiv, data = call_rate_indiv, FUN = sd)$rate

agg_rate_indiv <- agg_rate_indiv[agg_rate_indiv$indiv != "group", ]

mixed_ids <- unique(agg_rate_indiv$indiv[agg_rate_indiv$experiment == "mixed"])

mixed_l <- lapply(mixed_ids, function(x) {
  
  X <- agg_rate_indiv[agg_rate_indiv$indiv == x, ]
  X <- X[X$experiment == "mixed" | X$type == "solo", ]
  
  if (sum(X$experiment == "mixed") > 1)
    X <- X[X$group.size == max(X$group.size[X$experiment == "mixed"]) | X$type == "solo",]
  
  if (nrow(X) == 1)
    X <- X[vector(), ]
    
  return(X)
  })

mixed <- do.call(rbind, mixed_l)

mixed$experiment <- "mixed"

# regular
regular_ids <- unique(agg_rate_indiv$indiv[agg_rate_indiv$experiment == "regular"])

regular_l <- lapply(regular_ids, function(x) {
  
  X <- agg_rate_indiv[agg_rate_indiv$indiv == x, ]
  X <- X[X$experiment == "regular" | X$type == "solo", ]
  
  if (sum(X$experiment == "regular") > 1)
    X <- X[X$group.size == max(X$group.size[X$experiment == "regular"]) | X$type == "solo",]
  
  if (nrow(X) == 1)
    X <- X[vector(), ]
    
  return(X)
  })

regular <- do.call(rbind, regular_l)

regular$experiment <- "regular"

sub_agg <- rbind(mixed, regular)

sub_agg$experiment_f <- factor(ifelse(sub_agg$experiment == "regular", "Regular", "Artificial"), levels = c("Regular", "Artificial"))

sub_agg$type_f <- factor(ifelse(sub_agg$type == "solo", "Solo", "In group"), levels = c("Solo", "In group"))

# individual plot solo vs individual in group
ggplot(sub_agg, aes(x = type_f, y = rate, group = indiv, color = group.size)) + 
  geom_line() +
  geom_point(size = 2) +
  scale_color_viridis_c(alpha = 0.8) +
  labs(x = "Flight", y = "Call rate (calls/s)", color = "Group\nsize") + 
  facet_grid( ~ experiment_f) +
   theme(axis.text.x = element_text(angle = 45, hjust = 1))

Individual call rate in solo and group flight

call_rate_indiv <- call_rate_indiv[!is.na(call_rate_indiv$rate) & call_rate_indiv$rate > 0, ]

call_rate_indiv$experiment.type <- ifelse(call_rate_indiv$type == "solo", "solo", paste(call_rate_indiv$experiment, call_rate_indiv$type, sep = "-"))

call_rate_indiv$experiment.type <- gsub("mixed-indiv.group", "artificial.group", call_rate_indiv$experiment.type)
call_rate_indiv$experiment.type <- gsub("regular-indiv.group", "real.group", call_rate_indiv$experiment.type)

call_rate_indiv$experiment.type <- factor(call_rate_indiv$experiment.type, levels = c("solo", "real.group", "artificial.group")) 

# mean centering group size
call_rate_indiv$group.size <- call_rate_indiv$group.size - mean(call_rate_indiv$group.size)

 mod <- brm(
          formula = rate ~ experiment.type + experiment.type:group.size + (1 | indiv),
          iter = iter,
          thin = 1,
          data = call_rate_indiv,
          family = lognormal(),
          silent = 2,
          chains = chains,
          cores = chains, 
          prior = prior, 
          file = "./data/processed/individual_call_rate_solo_vs_group"
          )
html_summary(read.file = "./data/processed/individual_call_rate_solo_vs_group.rds", gsub.pattern = "experiment.type", gsub.replacement = "solo_vs_", remove.intercepts = TRUE, highlight = TRUE)

individual_call_rate_solo_vs_group

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20) rate ~ experiment.type + experiment.type:group.size + (1 | indiv) 10000 4 1 5000 0 0 26177.38 16125.3 852675893
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_solo_vs_real.group -1.297 -1.551 -1.044 1 34685.25 16697.46
b_solo_vs_artificial.group -1.520 -1.813 -1.223 1 30590.21 16125.30
b_solo_vs_solo:group.size -0.032 -0.139 0.074 1 26177.38 16509.06
b_solo_vs_real.group:group.size -0.139 -0.299 0.024 1 26201.68 16497.70
b_solo_vs_artificial.group:group.size -0.312 -0.532 -0.092 1 29927.08 16374.72

 

Takeaways

  • Individuals call at lower rates during group flight (both real and artificial) compared to solo flight
  • Rates also decreases in artificial groups as a function of group size

 

Individual call rate vs overall group call rate

call_rate_by_group <- readRDS("./data/processed/call_rate_by_group.RDS")

call_rate_by_group$experiment.type <- ifelse(call_rate_by_group$type == "solo", "solo", paste(call_rate_by_group$experiment, call_rate_by_group$type, sep = "-"))

call_rate_by_group$experiment.type <- gsub("mixed.group", "artificial.group", call_rate_by_group$experiment.type)
call_rate_by_group$experiment.type <- gsub("regular.group", "real.group", call_rate_by_group$experiment.type)

call_rate_by_group$experiment.type <- factor(call_rate_by_group$experiment.type, levels = c("solo", "real.group", "artificial.group")) 

# aggregate for plot
agg_rate_group <- aggregate(rate ~ experiment + type + group.size + group, data = call_rate_by_group, FUN = mean)

agg_rate_group$sd <- aggregate(rate ~ experiment + type + group.size + group, data = call_rate_by_group, FUN = sd)$rate

agg_rate_group <- agg_rate_group[agg_rate_group$type != "indiv.group", ]

mixed_ids <- unique(agg_rate_group$group[agg_rate_group$experiment == "mixed"])

mixed_l <- lapply(mixed_ids, function(x) {
  
  X <- agg_rate_group[agg_rate_group$group == x, ]
  X <- X[X$experiment == "mixed" | X$type == "solo", ]
  
  if (sum(X$experiment == "mixed") > 1)
    X <- X[X$group.size == max(X$group.size[X$experiment == "mixed"]) | X$type == "solo",]
  
  if (nrow(X) == 1)
    X <- X[vector(), ]
    
  return(X)
  })

mixed <- do.call(rbind, mixed_l)

mixed$experiment <- "mixed"

# regular
regular_ids <- unique(agg_rate_group$group[agg_rate_group$experiment == "regular"])

regular_l <- lapply(regular_ids, function(x) {
  
  X <- agg_rate_group[agg_rate_group$group == x, ]
  X <- X[X$experiment == "regular" | X$type == "solo", ]
  
  if (sum(X$experiment == "regular") > 1)
    X <- X[X$group.size == max(X$group.size[X$experiment == "regular"]) | X$type == "solo",]
  
  if (nrow(X) == 1)
    X <- X[vector(), ]
    
  return(X)
  })

regular <- do.call(rbind, regular_l)

regular$experiment <- "regular"

sub_agg <- rbind(mixed, regular)

sub_agg$experiment_f <- factor(ifelse(sub_agg$experiment == "regular", "Regular", "Artificial"), levels = c("Regular", "Artificial"))

sub_agg$type_f <- factor(ifelse(sub_agg$type == "solo", "Solo", "Overall group"), levels = c("Solo", "Overall group"))

# individual plot solo vs individual in group
ggplot(sub_agg, aes(x = type_f, y = rate, group = group, color = group.size)) + 
  geom_line() +
  geom_point(size = 2) +
  scale_color_viridis_c(alpha = 0.8) +
  labs(x = "Flight", y = "Call rate (calls/s)", color = "Group\nsize") + 
  facet_grid( ~ experiment_f) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# in the same scale than he one  above

# ggplot(sub_agg, aes(x = type_f, y = rate, group = group, color = group.size)) + 
#   geom_line() +
#   geom_point(size = 2) +
#   scale_color_viridis_c(alpha = 0.8) +
#   labs(x = "Flight", y = "Call rate (calls/s)", color = "Group\nsize") + 
#   facet_grid( ~ experiment_f) +
#   theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
#   ylim(c(0, 210))
call_rate_by_group <- readRDS("./data/processed/call_rate_by_group.RDS")

call_rate_by_group$experiment.type <- ifelse(call_rate_by_group$type == "solo", "solo", paste(call_rate_by_group$experiment, call_rate_by_group$type, sep = "-"))

call_rate_by_group$experiment.type <- gsub("mixed.group", "artificial.group", call_rate_by_group$experiment.type)
call_rate_by_group$experiment.type <- gsub("regular.group", "real.group", call_rate_by_group$experiment.type)

call_rate_by_group$experiment.type <- factor(call_rate_by_group$experiment.type, levels = c("solo", "real.group", "artificial.group")) 

# mean centering group size
call_rate_by_group$group.size <- call_rate_by_group$group.size - mean(call_rate_by_group$group.size)

 mod <- brm(
          formula = rate ~ experiment.type + experiment.type:group.size + (1 | group),
          iter = iter,
          thin = 1,
          data = call_rate_by_group,
          family = lognormal(),
          silent = 2,
          chains = chains,
          cores = chains, 
          prior = prior, 
          file = "./data/processed/individual_vs_group_call_rate"
          )
html_summary(read.file = "./data/processed/individual_vs_group_call_rate.rds", gsub.pattern = "experiment.type", gsub.replacement = "solo_vs_", remove.intercepts = TRUE, highlight = TRUE)

individual_vs_group_call_rate

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20) rate ~ experiment.type + experiment.type:group.size + (1 | group) 10000 4 1 5000 0 0 16437.37 14658.44 1799142864
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_solo_vs_real.group 0.292 -0.019 0.604 1 30238.76 14658.44
b_solo_vs_artificial.group 0.166 -0.159 0.487 1 22187.17 16136.54
b_solo_vs_solo:group.size -0.023 -0.109 0.068 1 16437.37 14988.15
b_solo_vs_real.group:group.size 0.023 -0.178 0.226 1 26569.51 14684.60
b_solo_vs_artificial.group:group.size 0.039 -0.224 0.300 1 21294.80 15810.40

 

Takeaways

  • No difference in overall call rate between solo and group flights (both real and artificial)

 

Gaps

Individual gaps in solo vs group flights

seltab <- attributes(clls)$check.results
seltab$year.audio <- clls$year.audio
seltab <- seltab[!is.na(seltab$sound.files), ]
seltab$type <- sapply(seltab$year.audio, function(x) metadat$Experimento[metadat$year.audio == x][1])
seltab$group <- sapply(seltab$year.audio, function(x) metadat$Grupo[metadat$year.audio == x][1])
seltab <- seltab[!is.na(seltab$group), ]
seltab$type <-ifelse(grepl("solo", seltab$type), "solo", "group")
seltab$indiv <- sapply(1:nrow(seltab), function(x) if (seltab$type[x] == "solo") metadat$Individuo[metadat$year.audio == seltab$year.audio[x]] else "group")
seltab$start <- seltab$orig.start
seltab$end <- seltab$orig.end

# get group and solo call counts per individual, call rate in call/min
gaps_indiv_list <- pblapply(1:nrow(diagnostics), function(x){

  # print(x)
  # get indivs
  sound.files <- acous_param_l[[diagnostics$group[x]]]$sound.files
  
    Y <- seltab[seltab$sound.files %in% sound.files, ]

  # add id for group flights
  Y$pred.id <- NA
  Y$pred.id[Y$type == "group"] <- sapply(Y$sound.files[Y$type == "group"], function(x) agg_pred$pred_indiv[agg_pred$sound.files == x])

  Y <- Y[order(Y$orig.sound.files, Y$start), ]
    
  gaps_indiv_list <- lapply(unique(Y$indiv), function(y){ 
# print(y)
    if (y != "group"){
      W <- Y[Y$indiv == y, ]
      W$sound.files <- W$orig.sound.files
      
      gaps.solo <- data.frame(group =  diagnostics$group[x], indiv = y, experiment =  diagnostics$experiment[x], type = "solo", gaps = if (nrow(W) > 0) na.omit(gaps(W, pb = FALSE)$gaps) else NA)
    
      Z <- Y[Y$pred.id == y & !is.na(Y$pred.id), ]
      Z$sound.files <- Z$orig.sound.files 
      group.gaps <- if (nrow(Z) > 0) na.omit(gaps(Z, pb = FALSE)$gaps) else NA
      gaps.group <- data.frame(group =  diagnostics$group[x], indiv = y, experiment =  diagnostics$experiment[x], type = "indiv.group", gaps =  if (length(group.gaps) > 0) group.gaps else NA)
       
    gaps_data <- rbind(gaps.group, gaps.solo)  
    
    } else {
      Q <- Y[Y$type == "group", ]
      Q$sound.files <- Q$orig.sound.files
    gaps_data <- data.frame(group =  diagnostics$group[x], indiv = y, experiment = diagnostics$experiment[x], type = "overall.group", gaps =  na.omit(gaps(Q, pb = FALSE)$gaps))
      
return(gaps_data)
    }
  })
  gaps_indiv <- do.call(rbind, gaps_indiv_list)
  gaps_indiv$group.size <- diagnostics$n_indiv[x]
  
  
    return(gaps_indiv)
    }
  )
  
gaps_indiv <- do.call(rbind, gaps_indiv_list)
  
# add sex
gaps_indiv$sex <- sapply(gaps_indiv$indiv, function(x) if(x != "group") na.exclude(caps$Sexo[caps$Murci == x])[1] else NA, USE.NAMES = FALSE)

gaps_indiv$sex <- ifelse(gaps_indiv$sex == "m", "Male", "Female")

# add age
gaps_indiv$age <- sapply(gaps_indiv$indiv, function(x) if(x != "group") na.exclude(caps$Edad[caps$Murci == x])[1] else NA, USE.NAMES = FALSE)

gaps_indiv$age <- ifelse(gaps_indiv$age == "sa", "Sub-adult", "Adult")

# reproductive stage
gaps_indiv$reprod.stg <- sapply(gaps_indiv$indiv, function(x) if(x != "group") na.exclude(caps$`Estado reproductivo`[caps$Murci == x])[1] else NA, USE.NAMES = FALSE)

gaps_indiv$reprod.stg[gaps_indiv$reprod.stg == "p?"] <- "Pregnant"
gaps_indiv$reprod.stg[gaps_indiv$reprod.stg == "ne"] <- "Inactive"

gaps_indiv$group.size <-sapply(gaps_indiv$group, function(x) diagnostics$n_indiv[diagnostics$group == x])

# order levels
gaps_indiv$experiment <- factor(gaps_indiv$experiment, levels = c("Regular", "Artificial"))

gaps_indiv$experiment.type <- ifelse(gaps_indiv$type == "solo", "solo", paste(gaps_indiv$experiment, gaps_indiv$type, sep = "-"))


# #find individuals with repeated solo data
# indivs_2_sets <- sapply(unique(gaps_indiv$indiv), function(x) length(unique(gaps_indiv$group[gaps_indiv$indiv == x & gaps_indiv$type == "solo"])))
# 
# # remove those
# indivs_2_sets <- names(indivs_2_sets)[indivs_2_sets > 1]
# indivs_2_sets <- indivs_2_sets[indivs_2_sets != "group"]
# 
# gaps_indiv_filter_list <- lapply(unique(gaps_indiv$indiv), function(x){
# 
#   Y <- gaps_indiv[gaps_indiv$indiv == x, ]
#   
#   if(x %in% indivs_2_sets){
#     Y <- Y[Y$group == unique(Y$group)[1] & Y$type == "solo" | Y$type == "indiv.group",]
#     
#   }
#   
#   return(Y)
#   
# })

# gaps_indiv_filter <- do.call(rbind, gaps_indiv_filter_list)
# 
# # all have 1
# all(sapply(unique(gaps_indiv_filter$indiv), function(x) length(unique(gaps_indiv_filter$group[gaps_indiv_filter$indiv == x & gaps_indiv_filter$type == "solo"]))) <= 1)

 saveRDS(gaps_indiv, "./data/processed/gaps_by_individual.RDS")
gaps_indiv <- readRDS("./data/processed/gaps_by_individual.RDS")


# aggregate for plot
agg_gaps <- aggregate(gaps ~ experiment + type + group.size, data = gaps_indiv, FUN = mean)

agg_gaps$sd <- aggregate(gaps ~ experiment + type + group.size, data = gaps_indiv, FUN = sd)$gaps

# ggplot(agg_gaps[agg_gaps$group.size < 6, ], aes(x = type, y = gaps)) + 
#   geom_point(color = viridis(10)[3], size = 2) +
#     geom_errorbar(aes(ymin = gaps - sd, ymax = gaps + sd), color = viridis(10)[3], width=.1, lwd = 1.1) +
#    geom_hline(yintercept = mean(gaps_indiv$gaps[gaps_indiv$type == "solo"], na.rm = TRUE), lty = 2, alpha = 0.5) + 
#   labs(x = "Flight", y = "Gap duration (s)") + 
#   facet_grid(group.size ~ experiment) +
#    theme(axis.text.x = element_text(angle = 45, hjust = 1))
# aggregate for plot
agg_gaps_indiv <- aggregate(gaps ~ experiment + type + group.size + indiv, data = gaps_indiv, FUN = mean)

agg_gaps_indiv$sd <- aggregate(gaps ~ experiment + type + group.size + indiv, data = gaps_indiv, FUN = sd)$gaps

agg_gaps_indiv <- agg_gaps_indiv[agg_gaps_indiv$indiv != "group", ]

mixed_ids <- unique(agg_gaps_indiv$indiv[agg_gaps_indiv$experiment == "mixed"])

mixed_l <- lapply(mixed_ids, function(x) {
  
  X <- agg_gaps_indiv[agg_gaps_indiv$indiv == x, ]
  X <- X[X$experiment == "mixed" | X$type == "solo", ]
  
  if (sum(X$experiment == "mixed") > 1)
    X <- X[X$group.size == max(X$group.size[X$experiment == "mixed"]) | X$type == "solo",]
  
  if (nrow(X) == 1)
    X <- X[vector(), ]
    
  return(X)
  })

mixed <- do.call(rbind, mixed_l)

mixed$experiment <- "mixed"

# regular
regular_ids <- unique(agg_gaps_indiv$indiv[agg_gaps_indiv$experiment == "regular"])

regular_l <- lapply(regular_ids, function(x) {
  
  X <- agg_gaps_indiv[agg_gaps_indiv$indiv == x, ]
  X <- X[X$experiment == "regular" | X$type == "solo", ]
  
  if (sum(X$experiment == "regular") > 1)
    X <- X[X$group.size == max(X$group.size[X$experiment == "regular"]) | X$type == "solo",]
  
  if (nrow(X) == 1)
    X <- X[vector(), ]
    
  return(X)
  })

regular <- do.call(rbind, regular_l)

regular$experiment <- "regular"

sub_agg <- rbind(mixed, regular)

sub_agg$experiment_f <- factor(ifelse(sub_agg$experiment == "regular", "Regular", "Artificial"), levels = c("Regular", "Artificial"))

sub_agg$type_f <- factor(ifelse(sub_agg$type == "solo", "Solo", "In group"), levels = c("Solo", "In group"))

# individual plot solo vs individual in group
ggplot(sub_agg, aes(x = type_f, y = gaps, group = indiv, color = group.size)) + 
  geom_line() +
  geom_point(size = 2) +
  scale_color_viridis_c(alpha = 0.8) +
  labs(x = "Flight", y = "Gap duration (s)", color = "Group\nsize") + 
  facet_grid( ~ experiment_f) +
   theme(axis.text.x = element_text(angle = 45, hjust = 1))

sub_gaps_indiv <- gaps_indiv[gaps_indiv$type != "overall.group", ]

sub_gaps_indiv$experiment.type <- gsub("-indiv.group", "", sub_gaps_indiv$experiment.type)

sub_gaps_indiv$experiment.type <- gsub("mixed", "artificial.group", sub_gaps_indiv$experiment.type)
sub_gaps_indiv$experiment.type <- gsub("regular", "real.group", sub_gaps_indiv$experiment.type)

sub_gaps_indiv$experiment.type <- factor(sub_gaps_indiv$experiment.type, levels = c("solo", "real.group", "artificial.group")) 

# mean centering group size
sub_gaps_indiv$group.size <- sub_gaps_indiv$group.size - mean(sub_gaps_indiv$group.size)

 mod <- brm(
          formula =  gaps ~ experiment.type + experiment.type:group.size + (1 | indiv),
          iter = iter,
          thin = 1,
          data = sub_gaps_indiv,
          family = lognormal(),
          silent = 2,
          chains = chains,
          cores = chains, 
          prior = prior, 
          file = "./data/processed/individual_gaps_solo_vs_group"
          )
html_summary(read.file = "./data/processed/individual_gaps_solo_vs_group.rds", gsub.pattern = "experiment.type", gsub.replacement = "solo_vs_", remove.intercepts = TRUE, highlight = TRUE)

individual_gaps_solo_vs_group

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20) gaps ~ experiment.type + experiment.type:group.size + (1 | indiv) 10000 4 1 5000 0 0 4768.63 9070.179 1328396827
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_solo_vs_real.group 0.444 0.376 0.510 1 16603.356 15045.895
b_solo_vs_artificial.group 0.650 0.576 0.723 1 10885.827 13391.265
b_solo_vs_solo:group.size -0.035 -0.072 0.001 1.001 4768.630 9070.179
b_solo_vs_real.group:group.size -0.041 -0.097 0.015 1 8557.013 12716.294
b_solo_vs_artificial.group:group.size 0.174 0.123 0.225 1 11901.013 14209.919

# contrasts
cat("Compare artificial vs real group gaps:\n")

Compare artificial vs real group gaps:

# contrasts
contrasts(read.file = "./data/processed/individual_gaps_solo_vs_group.rds", predictor = "experiment.type", n.posterior = 2000, level.sep = " VS ",  html.table = TRUE, plot = TRUE, highlight = TRUE)
Hypothesis Estimate Est.Error l-95% CI u-95% CI
1 real.group VS solo 0.444 0.034 0.376 0.510
2 real.group VS artificial.group -0.206 0.046 -0.296 -0.117
3 solo VS artificial.group -0.650 0.038 -0.576 -0.723

sub_gaps_indiv <- gaps_indiv[gaps_indiv$type != "overall.group", ]

sub_gaps_indiv$experiment.type <- gsub("-indiv.group", "", sub_gaps_indiv$experiment.type)

sub_gaps_indiv$experiment.type <- gsub("mixed", "artificial.group", sub_gaps_indiv$experiment.type)
sub_gaps_indiv$experiment.type <- gsub("regular", "real.group", sub_gaps_indiv$experiment.type)

sub_gaps_indiv$experiment.type <- factor(sub_gaps_indiv$experiment.type, levels = c("solo", "real.group", "artificial.group")) 

# mean centering group size
sub_gaps_indiv$group.size <- sub_gaps_indiv$group.size - mean(sub_gaps_indiv$group.size)


 mod <- brm(
          formula = gaps ~ experiment.type + experiment.type:group.size + (experiment.type | indiv),
          iter = iter,
          thin = 1,
          data = sub_gaps_indiv,
          family = lognormal(),
          silent = 2,
          chains = chains,
          cores = chains, 
          prior = prior, 
          file = "./data/processed/individual_gaps_solo_vs_group_random_slope"
          )
html_summary(read.file = "./data/processed/individual_gaps_solo_vs_group_random_slope.rds", gsub.pattern = "experiment.type", gsub.replacement = "solo_vs_", remove.intercepts = TRUE, highlight = TRUE)

# 
# # contrasts
cat("Compare artificial vs real group gaps:\n")
# contrasts
contrasts(read.file = "../data/processed/individual_gaps_solo_vs_group_random_slope.rds", predictor = "experiment.type", n.posterior = 2000, level.sep = " VS ",  html.table = TRUE, plot = TRUE, highlight = TRUE)

 

Takeaways

  • Gaps in group flights (both mixed and regular groups) increases compare to solo flights but the magnitude of the response varies individually
  • Group size also affects the magnitude of change in group flights

 

Individual gaps vs overall group gaps

# aggregate for plot
agg_gaps_group <- aggregate(gaps ~ experiment + type + group.size + group, data = gaps_indiv, FUN = mean)

agg_gaps_group$sd <- aggregate(gaps ~ experiment + type + group.size + group, data = gaps_indiv, FUN = sd)$gaps

agg_gaps_group <- agg_gaps_group[agg_gaps_group$type != "indiv.group", ]

mixed_ids <- unique(agg_gaps_group$group[agg_gaps_group$experiment == "mixed"])

mixed_l <- lapply(mixed_ids, function(x) {
  
  X <- agg_gaps_group[agg_gaps_group$group == x, ]
  X <- X[X$experiment == "mixed" | X$type == "solo", ]
  
  if (sum(X$experiment == "mixed") > 1)
    X <- X[X$group.size == max(X$group.size[X$experiment == "mixed"]) | X$type == "solo",]
  
  if (nrow(X) == 1)
    X <- X[vector(), ]
    
  return(X)
  })

mixed <- do.call(rbind, mixed_l)

mixed$experiment <- "mixed"

# regular
regular_ids <- unique(agg_gaps_group$group[agg_gaps_group$experiment == "regular"])

regular_l <- lapply(regular_ids, function(x) {
  
  X <- agg_gaps_group[agg_gaps_group$group == x, ]
  X <- X[X$experiment == "regular" | X$type == "solo", ]
  
  if (sum(X$experiment == "regular") > 1)
    X <- X[X$group.size == max(X$group.size[X$experiment == "regular"]) | X$type == "solo",]
  
  if (nrow(X) == 1)
    X <- X[vector(), ]
    
  return(X)
  })

regular <- do.call(rbind, regular_l)

regular$experiment <- "regular"

sub_agg <- rbind(mixed, regular)

sub_agg$experiment_f <- factor(ifelse(sub_agg$experiment == "regular", "Regular", "Artificial"), levels = c("Regular", "Artificial"))

sub_agg$type_f <- factor(ifelse(sub_agg$type == "solo", "Solo", "Overall group"), levels = c("Solo", "Overall group"))

# individual plot solo vs individual in group
ggplot(sub_agg, aes(x = type_f, y = gaps, group = group, color = group.size)) + 
  geom_line() +
  geom_point(size = 2) +
  scale_color_viridis_c(alpha = 0.8) +
  labs(x = "Flight", y = "Gap duration (s)", color = "Group\nsize") + 
  facet_grid( ~ experiment_f) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# in the same scale than he one  above

ggplot(sub_agg, aes(x = type_f, y = gaps, group = group, color = group.size)) + 
  geom_line() +
  geom_point(size = 2) +
  scale_color_viridis_c(alpha = 0.8) +
  labs(x = "Flight", y = "Gap duration (s)", color = "Group\nsize") + 
  facet_grid( ~ experiment_f) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
  ylim(c(0, 210))

sub_gaps_indiv <- gaps_indiv[gaps_indiv$type != "indiv.group", ]

sub_gaps_indiv$experiment.type <- gsub("-overall.group", "", sub_gaps_indiv$experiment.type)

sub_gaps_indiv$experiment.type <- gsub("mixed", "artificial.group", sub_gaps_indiv$experiment.type)
sub_gaps_indiv$experiment.type <- gsub("regular", "real.group", sub_gaps_indiv$experiment.type)

sub_gaps_indiv$experiment.type <- factor(sub_gaps_indiv$experiment.type, levels = c("solo", "real.group", "artificial.group")) 

 mod <- brm(
          formula =  gaps ~ experiment.type + experiment.type:group.size + (1 | group),
          iter = iter,
          thin = 1,
          data = sub_gaps_indiv,
          family = lognormal(),
          silent = 2,
          chains = chains,
          cores = chains, 
          prior = prior, 
          file = "./data/processed/group_gaps_solo_vs_group"
 )
html_summary(read.file = "./data/processed/group_gaps_solo_vs_group.rds", gsub.pattern = "experiment.type", gsub.replacement = "solo_vs_", remove.intercepts = TRUE, highlight = TRUE)

group_gaps_solo_vs_group

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20) gaps ~ experiment.type + experiment.type:group.size + (1 | group) 10000 4 1 5000 0 0 1395.009 2988.005 790846726
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_solo_vs_real.group 0.141 -0.058 0.337 1 14050.802 14283.191
b_solo_vs_artificial.group 0.346 0.086 0.606 1 8320.936 12260.152
b_solo_vs_solo:group.size -0.014 -0.108 0.081 1.001 1395.009 2988.005
b_solo_vs_real.group:group.size -0.187 -0.286 -0.086 1.001 1510.895 3321.743
b_solo_vs_artificial.group:group.size -0.169 -0.274 -0.063 1.001 1680.001 3767.807

# contrasts
cat("Compare artificial vs real group gaps:\n")

Compare artificial vs real group gaps:

# contrasts
contrasts(read.file = "./data/processed/group_gaps_solo_vs_group.rds", predictor = "experiment.type", n.posterior = 2000, level.sep = " VS ",  html.table = TRUE, plot = TRUE, highlight = TRUE)
Hypothesis Estimate Est.Error l-95% CI u-95% CI
1 solo VS real.group -0.141 0.101 0.058 -0.337
2 solo VS artificial.group -0.346 0.132 -0.086 -0.606
3 real.group VS artificial.group -0.205 0.165 -0.526 0.116

Time partitioning during group flight

eltab <- attributes(clls)$check.results
seltab$year.audio <- clls$year.audio
seltab <- seltab[!is.na(seltab$sound.files), ]
seltab$type <- sapply(seltab$year.audio, function(x) metadat$Experimento[metadat$year.audio == x][1])
seltab$group <- sapply(seltab$year.audio, function(x) metadat$Grupo[metadat$year.audio == x][1])
seltab <- seltab[!is.na(seltab$group), ]
# seltab$type <-ifelse(grepl("solo", seltab$type), "solo", "group")
seltab$indiv <- sapply(1:nrow(seltab), function(x) if (seltab$type[x] == "solo") metadat$Individuo[metadat$year.audio == seltab$year.audio[x]] else "group")
seltab$start <- seltab$orig.start
seltab$end <- seltab$orig.end

# get group and solo call counts per individual, call rate in call/min
pred_indiv_list <- pblapply(1:nrow(diagnostics), function(x){

  # print(x)
  # get indivs
  sound.files <- acous_param_l[[diagnostics$group[x]]]$sound.files
  
    Y <- seltab[seltab$sound.files %in% sound.files, ]
    Y <- Y[grep("grup", Y$type), ]
  # add id for group flights
  Y$pred.id <- NA
  Y$pred.id <- sapply(Y$sound.files, function(x){ 
    w <- agg_pred$pred_indiv[agg_pred$sound.files == x]
    if (length(w) == 0) w <- NA
    return(w)
    })

  Y <- Y[order(Y$orig.sound.files, Y$start), ]
  Y$start <- Y$orig.start
  Y$end <- Y$orig.end
  Y$selec <- Y$orig.selec
  Y$sing.event <- Y$orig.sound.files
  Y$indiv <- Y$pred.id
  return(Y[, c("sing.event", "indiv", "selec", "start", "end", "group", "type")])
    }
  )
  
pred_indiv_df <- do.call(rbind, pred_indiv_list)
  
sum(is.na(pred_indiv_df$indiv)) / nrow(pred_indiv_df)

unique(pred_indiv_df$x[is.na(pred_indiv_df$indiv)])


pred_indiv_df <- pred_indiv_df[!is.na(pred_indiv_df$indiv), ]

# add sex
pred_indiv_df$sex <- sapply(pred_indiv_df$indiv, function(x) if(x != "group") na.exclude(caps$Sexo[caps$Murci == x])[1] else NA, USE.NAMES = FALSE)

pred_indiv_df$sex <- ifelse(pred_indiv_df$sex == "m", "Male", "Female")

# add age
pred_indiv_df$age <- sapply(pred_indiv_df$indiv, function(x) if(x != "group") na.exclude(caps$Edad[caps$Murci == x])[1] else NA, USE.NAMES = FALSE)

pred_indiv_df$age <- ifelse(pred_indiv_df$age == "sa", "Sub-adult", "Adult")

# reproductive stage
pred_indiv_df$reprod.stg <- sapply(pred_indiv_df$indiv, function(x) if(x != "group") na.exclude(caps$`Estado reproductivo`[caps$Murci == x])[1] else NA, USE.NAMES = FALSE)

pred_indiv_df$reprod.stg[pred_indiv_df$reprod.stg == "p?"] <- "Pregnant"
pred_indiv_df$reprod.stg[pred_indiv_df$reprod.stg == "ne"] <- "Inactive"

pred_indiv_df$group.size <-sapply(pred_indiv_df$group, function(x) diagnostics$n_indiv[diagnostics$group == x])

# order levels
pred_indiv_df$experiment <- ifelse(grepl("mixto", pred_indiv_df$type), "Artificial", "Regular")

pred_indiv_df$experiment <- factor(pred_indiv_df$experiment, levels = c("Regular", "Artificial"))

source("~/Dropbox/R_package_testing/warbleR/R/test_coordination.R")
source("~/Dropbox/R_package_testing/warbleR/R/internal_functions.R")

coordination_results <- test_coordination(X = pred_indiv_df, iterations = 10000, less.than.chance = TRUE, ovlp.method = "time.closest", parallel = 20, cutoff = 5)

coordination_results$n.calls <- sapply(coordination_results$sing.event, function(x) sum(pred_indiv_df$sing.event == x))

coordination_results$n.calls <- sapply(coordination_results$sing.event, function(x) sum(pred_indiv_df$sing.event == x))

coordination_results$experiment <- sapply(coordination_results$sing.event, function(x) pred_indiv_df$experiment[pred_indiv_df$sing.event == x][1])

coordination_results$group.size <- sapply(coordination_results$sing.event, function(x) length(unique(pred_indiv_df$indiv[pred_indiv_df$sing.event == x])))

coordination_results <- coordination_results[!is.na(coordination_results$p.value), ]

write.csv(coordination_results, "./data/processed/vocal_coordination_test_results.csv", row.names = FALSE)
coordination_results <- read.csv("./data/processed/vocal_coordination_test_results.csv")

coordination_results <- coordination_results[!is.na(coordination_results$coor.score),]


hist(coordination_results$p.value, breaks = 20)
hist(coordination_results$coor.score, breaks = 20)



coordination_results$coor.score_cr <- residuals(lm(coor.score ~ group.size, coordination_results))


ggplot(data = coordination_results, mapping = aes(x = group.size, y = p.value)) +
  # geom_hline(yintercept = 0, col = "red") +
  geom_point() + 
  labs(x = "Group size", y = "p value") +
  theme_classic()

ggplot(data = coordination_results, mapping = aes(x = group.size, y = coor.score)) +
  # geom_hline(yintercept = 0, col = "red") +
  geom_point() + 
  labs(x = "Group size", y = "Coordination score") +
  theme_classic()

ggplot(data = coordination_results, mapping = aes(x = group.size, y = coor.score)) +
  # geom_hline(yintercept = 0, col = "red") +
  geom_point() + 
  labs(x = "Group size", y = "Coordination score") +
  theme_classic()

ggplot(data = coordination_results, mapping = aes(x = group.size, y = coor.score_cr)) +
  # geom_hline(yintercept = 0, col = "red") +
  geom_point() + 
  labs(x = "Group size", y = "Coordination score") +
  theme_classic()

ggplot(data = coordination_results, mapping = aes(x = experiment, y = coor.score_cr)) +
  geom_hline(yintercept = 0, col = "red") +
  geom_boxplot() + 
  labs(x = "Group", y = "Coordination score") +
  theme_classic()
  
metap::sumz(1 - coordination_results$p.value,  weights = coordination_results$group.size)

metap::sumz(1 - coordination_results$p.value[coordination_results$experiment == "Regular"], weights = coordination_results$group.size[coordination_results$experiment == "Regular"])

metap::sumz(1 - coordination_results$p.value[coordination_results$experiment == "Artificial"], weights = coordination_results$group.size[coordination_results$experiment == "Artificial"])

mod <- brm(
          formula = coor.score_cr ~experiment,
          iter = iter,
          thin = 1,
          data = coordination_results,
          family = gaussian(),
          silent = 2,
          chains = chains,
          cores = chains,
          prior = priors, 
          file = "./data/processed/regression_coordination_score_by_group_type")
)  

 

Takeaways

  • Gaps in artificial groups are longer than in solo flights (but not for real groups)
  • Group size positively affects gap duration in both types of group flights

 


Session information

## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
## 
## locale:
##  [1] LC_CTYPE=pt_BR.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=pt_BR.UTF-8    
##  [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=pt_BR.UTF-8   
##  [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] metap_1.8           multtest_2.54.0     Biobase_2.54.0     
##  [4] BiocGenerics_0.40.0 purrr_1.0.0         tidybayes_3.0.1    
##  [7] ggridges_0.5.4      posterior_1.3.1     ggrepel_0.9.1      
## [10] brms_2.18.0         Rcpp_1.0.9          soundgen_2.1.0     
## [13] shinyBS_0.61        Sim.DiffProc_4.8    rlang_1.0.6        
## [16] kableExtra_1.3.4    DT_0.26             viridis_0.6.2      
## [19] viridisLite_0.4.1   pbapply_1.6-0       e1071_1.7-7        
## [22] caret_6.0-88        ggplot2_3.4.0       lattice_0.20-44    
## [25] ranger_0.13.1       readxl_1.3.1        warbleR_1.1.28     
## [28] NatureSounds_1.0.4  knitr_1.42          seewave_2.2.0      
## [31] tuneR_1.4.1         devtools_2.4.5      usethis_2.1.6      
## [34] lubridate_1.7.10    data.table_1.14.0  
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.2           tidyselect_1.2.0     lme4_1.1-27.1       
##   [4] fftw_1.0-7           htmlwidgets_1.5.4    grid_4.1.0          
##   [7] pROC_1.17.0.1        munsell_0.5.0        mutoss_0.1-12       
##  [10] codetools_0.2-18     miniUI_0.1.1.1       withr_2.5.0         
##  [13] Brobdingnag_1.2-9    colorspace_2.0-3     highr_0.10          
##  [16] rstudioapi_0.14      stats4_4.1.0         dtw_1.23-1          
##  [19] bayesplot_1.9.0      labeling_0.4.2       Rdpack_2.4          
##  [22] emmeans_1.8.1-1      rstan_2.21.7         mnormt_2.0.2        
##  [25] farver_2.1.1         bridgesampling_1.1-2 coda_0.19-4         
##  [28] vctrs_0.5.2          generics_0.1.3       TH.data_1.1-0       
##  [31] ipred_0.9-11         xfun_0.36            R6_2.5.1            
##  [34] markdown_1.3         gamm4_0.2-6          projpred_2.0.2      
##  [37] bitops_1.0-7         cachem_1.0.6         assertthat_0.2.1    
##  [40] promises_1.2.0.1     scales_1.2.1         multcomp_1.4-17     
##  [43] nnet_7.3-16          gtable_0.3.1         processx_3.8.0      
##  [46] sandwich_3.0-1       timeDate_3043.102    systemfonts_1.0.4   
##  [49] splines_4.1.0        ModelMetrics_1.2.2.2 checkmate_2.1.0     
##  [52] inline_0.3.19        yaml_2.3.7           reshape2_1.4.4      
##  [55] abind_1.4-5          threejs_0.3.3        crosstalk_1.2.0     
##  [58] backports_1.4.1      httpuv_1.6.6         tensorA_0.36.2      
##  [61] tools_4.1.0          lava_1.6.9           ellipsis_0.3.2      
##  [64] jquerylib_0.1.4      proxy_0.4-27         sessioninfo_1.2.2   
##  [67] TFisher_0.2.0        plyr_1.8.7           base64enc_0.1-3     
##  [70] RCurl_1.98-1.9       ps_1.7.2             prettyunits_1.1.1   
##  [73] rpart_4.1-15         cowplot_1.1.1        urlchecker_1.0.1    
##  [76] zoo_1.8-11           fs_1.6.0             magrittr_2.0.3      
##  [79] ggdist_3.2.0         colourpicker_1.2.0   tmvnsim_1.0-2       
##  [82] mvtnorm_1.1-3        matrixStats_0.62.0   pkgload_1.3.2       
##  [85] arrayhelpers_1.1-0   shinyjs_2.1.0        mime_0.12           
##  [88] evaluate_0.20        xtable_1.8-4         shinystan_2.6.0     
##  [91] gridExtra_2.3        rstantools_2.2.0     compiler_4.1.0      
##  [94] tibble_3.1.8         crayon_1.5.2         minqa_1.2.4         
##  [97] StanHeaders_2.21.0-7 htmltools_0.5.4      mgcv_1.8-36         
## [100] later_1.3.0          tidyr_1.1.3          RcppParallel_5.1.5  
## [103] DBI_1.1.1            MASS_7.3-54          boot_1.3-28         
## [106] Matrix_1.5-1         cli_3.6.0            rbibutils_2.2.9     
## [109] qqconf_1.3.1         parallel_4.1.0       gower_0.2.2         
## [112] igraph_1.3.5         pkgconfig_2.0.3      sn_2.1.0            
## [115] numDeriv_2016.8-1.1  signal_0.7-7         recipes_0.1.16      
## [118] svUnit_1.0.6         xml2_1.3.3           foreach_1.5.1       
## [121] dygraphs_1.1.1.6     svglite_2.1.0        bslib_0.4.2         
## [124] webshot_0.5.4        estimability_1.4.1   prodlim_2019.11.13  
## [127] rvest_1.0.3          stringr_1.5.0        distributional_0.3.1
## [130] callr_3.7.3          digest_0.6.31        rmarkdown_2.20      
## [133] cellranger_1.1.0     Deriv_4.1.3          shiny_1.7.3         
## [136] gtools_3.9.3         rjson_0.2.21         nloptr_1.2.2.2      
## [139] lifecycle_1.0.3      nlme_3.1-152         jsonlite_1.8.4      
## [142] fansi_1.0.3          pillar_1.8.1         loo_2.4.1.9000      
## [145] plotrix_3.8-1        fastmap_1.1.0        httr_1.4.4          
## [148] pkgbuild_1.4.0       survival_3.2-11      glue_1.6.2          
## [151] xts_0.12.2           remotes_2.4.2        shinythemes_1.2.0   
## [154] iterators_1.0.13     class_7.3-19         stringi_1.7.12      
## [157] sass_0.4.5           profvis_0.3.7        memoise_2.0.1       
## [160] mathjaxr_1.6-0       dplyr_1.0.10